Progress Made

Group Roles


Research Question 1

Research Question 2

Load Libraries

library(tidyr)
library(rjags)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(MacBayes)
library(choroplethrZip)
library(rgdal)

Map cleanup

cleanup <- theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = 'white', colour = 'white'),
        axis.line = element_line(colour = "white"),
        axis.ticks=element_blank(), axis.text.x=element_blank(),
        axis.text.y=element_blank())

Clean data

#Load data
Crime = read.csv('Crime_sample_2010_to_2017.csv')

Crime = Crime %>%
  mutate(Date = as.POSIXct(strptime(Date, "%Y-%m-%d %H:%M:%S"))) %>%
  filter(!is.na(Latitude))
## Warning in strptime(Date, "%Y-%m-%d %H:%M:%S"): unknown timezone 'zone/tz/
## 2017c.1.0/zoneinfo/America/Chicago'

Project observations to ZIP codes boundaries

#load shapefile
zip_shp <- readOGR(dsn="Boundaries - ZIP Codes", layer="geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a")
## OGR data source with driver: ESRI Shapefile 
## Source: "Boundaries - ZIP Codes", layer: "geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a"
## with 61 features
## It has 4 fields
temp <- zip_shp
zip_df <- fortify(temp, region = "zip")

#needs to reassign CRS for shapefile
new_CRS <- CRS("+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
zip_shp <- sp::spTransform(zip_shp, new_CRS)

#summarize crime count per year
years <- unique(Crime$Year)
Hotspot_Crime <- data.frame()

for(i in 1:length(years)){
  year <- years[i]
  Current_Year <- Crime %>% filter(Year == year)
  
  locations <- with(Current_Year, as.data.frame(cbind(Longitude, Latitude)))
  coordinates(locations) <- ~Longitude+Latitude
  proj4string(locations) <- CRS("+init=epsg:4326")
      
  #prepping the data: project points to polygons
  proj4string(zip_shp) <- new_CRS
  proj4string(locations) <- new_CRS
  by_zip <- over(locations, zip_shp)
  
  by_zip <- by_zip %>%
    group_by(zip) %>%
    dplyr::summarise(total=n()) %>%
    filter(!is.na(zip)) %>%
    mutate(id = as.character(zip))
  
  total_map <- left_join(zip_df, by_zip)
  total_map <- total_map %>% mutate(Year = year)
  
  Hotspot_Crime <- rbind(Hotspot_Crime, total_map)
}

Add demographic information to Model Data

data("df_zip_demographics")
Demographics <- df_zip_demographics %>% mutate(zip = region)
Hotspot_Crime <- inner_join(Hotspot_Crime, Demographics, by="zip") 
## Warning: Column `zip` joining factor and character vector, coercing into
## character vector
Model_Data <- Hotspot_Crime %>%
  select(-c(lat, long, order, piece, group, hole)) %>% distinct

Get populations by zipcode and year

Crime_count <- Model_Data %>%
  select(Year, zip, total) %>%
  spread(key=Year, value=total)
rownames(Crime_count) <- Crime_count$zip 
Crime_count <- Crime_count %>% select(-zip)
Crime_count[is.na(Crime_count)] <- 0
head(Crime_count)
##       2010 2011 2012 2013 2014 2015 2016 2017
## 60601    8    6    5    8    7    6    7   11
## 60602    4    8    1    4    5   10    8    9
## 60603    1    5    3    2    1    2    9    4
## 60604    1    6    3    1    4    3    1    3
## 60605   10    7   13    6    9   10   10    8
## 60606    2    6    4    3    1    4    3    2
Population <- Model_Data %>%
  select(total_population, zip, Year) %>%
  spread(key=Year, value=total_population) %>%
  select('2010', zip)
rownames(Population) <- Population$zip
Population <- Population %>% select(-zip)
Population[is.na(Population)] <- 0
head(Population)
##        2010
## 60601  9732
## 60602  1463
## 60603   882
## 60604   415
## 60605 24972
## 60606  2789
Population.vector <- Population$`2010`
time <- sort(unique(Model_Data$Year))
percent_nonwhite <- 100 - Model_Data$percent_white

Simulation

Research Question 1

Plot of \(log(y_{ij}/n_{i})\) vs year

Crime_count_log_calc <- log(Crime_count/Population.vector)

Crime_count_log <- Crime_count_log_calc %>%
  rownames_to_column() %>%
  gather(key=year, value=crime_count_log, -rowname) %>%
  mutate(zip=rowname) %>%
  select(zip, year, crime_count_log) %>%
  mutate(year=as.numeric(year))
head(Crime_count_log)
##     zip year crime_count_log
## 1 60601 2010       -7.103733
## 2 60602 2010       -5.901950
## 3 60603 2010       -6.782192
## 4 60604 2010       -6.028279
## 5 60605 2010       -7.822925
## 6 60606 2010       -7.240291
ggplot(Crime_count_log, aes(x=year, y=crime_count_log, colour=zip)) + geom_point() + geom_line() 

Plot of \(log(y_{ij}/n_{i})\) vs per capita income

Model_Data_1 <- Model_Data %>%
  filter(Year == 2010) %>%
  mutate(percent_nonwhite = 100 - percent_white) %>%
  select(zip, per_capita_income, percent_nonwhite) %>%
  right_join(Crime_count_log, by="zip") 
head(Model_Data_1)
##     zip per_capita_income percent_nonwhite year crime_count_log
## 1 60601             89136               37 2010       -7.103733
## 2 60602             66748               38 2010       -5.901950
## 3 60603             73239               39 2010       -6.782192
## 4 60604            135807               20 2010       -6.028279
## 5 60605             58338               42 2010       -7.822925
## 6 60606             83005               25 2010       -7.240291
ggplot(Model_Data_1, aes(x=log(per_capita_income), y=crime_count_log, colour=percent_nonwhite)) + geom_point() + labs(x="Log of Per Capita Income", y="Log of Crime Rate", title="Relationship between Per Capita Income and Log of Crime Rate")

Plot of \(log(y_{ij}/n_{i})\) vs nonwhite percentage

ggplot(Model_Data_1, aes(x=log(percent_nonwhite), y=crime_count_log)) + geom_point() + labs(x="Log of Nonwhite Percentage", y="Log of Crime Rate", title="Relationship between Nonwhite Percentage and Log of Crime Rate")

Model 0

Data

\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \beta * time_{j}\]

Priors

\[\beta \sim N(0, 1000^2)\]

  • i: zipcodes in Chicago (58 total in our data)

  • j: years (from 2010 to 2017)

  • y[i, j]: observed count of violent crimes in Chicago by year and zipcode

  • p[i, j]: probability of being a victim in location i and time j

  • n[i]: the population in zipcode i

  • \(\beta\): time trend

  • time[j]: the time in years

#Data
Crime_count <- Model_Data %>%
  group_by(Year) %>%
  dplyr::summarise(total.stop = sum(total))
head(Crime_count)
## # A tibble: 6 x 2
##    Year total.stop
##   <int>      <int>
## 1  2010       1474
## 2  2011       1448
## 3  2012       1340
## 4  2013       1258
## 5  2014       1090
## 6  2015       1117
#Specify the model
crime_model_0 <- "model{
    #Data
    for(j in 1:8){
      y[j] ~ dbin(p[j], n)
      logit(p[j]) = beta*time[j]
    }
    
    #Priors
    beta ~ dnorm(0, 1/100^2)

}"

#set up an algorithm to simulate the posterior by 
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_0 <- jags.model(textConnection(crime_model_0), data=list(y=Crime_count$total.stop, n=sum(Population.vector), time=Crime_count$Year),inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 8
##    Unobserved stochastic nodes: 1
##    Total graph size: 40
## 
## Initializing model
#simulate a sample from the posterior 
#note that we specify both mu and tau variables
crime_sim_0 <- coda.samples(crime_jags_0, variable.names=c("p", "beta"), n.iter=50000)

#store the samples in a data frame:    
crime_data_0 <- data.frame(crime_sim_0[[1]])

Running Mean Plot

running_mean_plot(crime_data_0$beta)

Parameter Inferences

ggplot(crime_data_0, aes(x=beta)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Effect of Time", x="Time Effect", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_0, aes(x=p.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2010", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_0, aes(x=p.2.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2011", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_0, aes(x=p.3.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2012", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_0, aes(x=p.4.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2013", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_0, aes(x=p.5.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2014", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_0, aes(x=p.6.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2015", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_0, aes(x=p.7.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2016", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_0, aes(x=p.8.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2017", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Confidence Intervals and Posterior Means of Crime Model 0 Parameters

crime_model_0_CI <- matrix(nrow=9, ncol=2)

for (i in 1:9) {
      crime_model_0_CI[i,] <- quantile(crime_data_0[[i]], c(0.025, 0.975))
}

crime_model_0_mean <- matrix(nrow=9, ncol=1)
for (i in 1:9) {
      crime_model_0_mean[i, ] <- mean(crime_data_0[[i]])
}

crime_model_0_table <- cbind(crime_model_0_CI, crime_model_0_mean)
rownames(crime_model_0_table) <- names(crime_data_0)
colnames(crime_model_0_table) <- c("2.5%", "97.5%", "Mean")
as.table(crime_model_0_table)
##               2.5%         97.5%          Mean
## beta -0.0038479706 -0.0038283390 -0.0038380756
## p.1.  0.0004373144  0.0004549075  0.0004461179
## p.2.  0.0004356356  0.0004531701  0.0004444097
## p.3.  0.0004339632  0.0004514393  0.0004427081
## p.4.  0.0004322972  0.0004497151  0.0004410130
## p.5.  0.0004306377  0.0004479975  0.0004393244
## p.6.  0.0004289845  0.0004462865  0.0004376422
## p.7.  0.0004273377  0.0004445820  0.0004359665
## p.8.  0.0004256971  0.0004428839  0.0004342971

Model 1

Notation

  • i: zipcodes in Chicago (58 total in our data)

  • j: years (from 2010 to 2017)

  • y[i, j]: observed count of violent crimes in Chicago by year and zipcode

  • p[i, j]: probability of being a victim in location i and time j

  • n[i, j]: the population in zipcode i in time j.

  • time[j]: the time in years

  • \(\delta\)[i]: the spatial temporal interaction term

Specify the Model

We are modeling the count of crime incidents by year and zipcodes using a binomial distribution, in which the parameters are the probability of being a victim and the total population in location i and year j.

We modeled the probability of being a victim using a logit function of the grand crime rate of Chicago and the inclusion of several random effects. It is also a function of time and its interaction with location.

Data

\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \delta_{i} * time_{j}\]

Priors

\[\delta_{i} \sim N(0, prec.delta)\]

Hyperpriors

\[prec.delta \sim Gamma(0.5, 0.0005)\]

#Data
Crime_count <- Model_Data %>%
  select(Year, zip, total) %>%
  spread(key=Year, value=total)
rownames(Crime_count) <- Crime_count$zip 
Crime_count <- Crime_count %>% select(-zip)
Crime_count[is.na(Crime_count)] <- 0
head(Crime_count)
##       2010 2011 2012 2013 2014 2015 2016 2017
## 60601    8    6    5    8    7    6    7   11
## 60602    4    8    1    4    5   10    8    9
## 60603    1    5    3    2    1    2    9    4
## 60604    1    6    3    1    4    3    1    3
## 60605   10    7   13    6    9   10   10    8
## 60606    2    6    4    3    1    4    3    2
#Specify the model
crime_model_1 <- "model{
    #Data
    for(i in 1:58) {
      for(j in 1:8){
        y[i,j] ~ dbin(p[i,j], n[i])
        logit(p[i,j]) = delta[i]*time[j]
      }
    }
    
    #Priors
    for(i in 1:58){
      delta[i] ~ dnorm(0, prec.delta)
    }

    #Hyperpriors
    prec.delta ~ dgamma(0.5, 0.0005)
    
}"

#set up an algorithm to simulate the posterior by 
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_1 <- jags.model(textConnection(crime_model_1), data=list(y=Crime_count, n=Population.vector, time=time),
                    inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 464
##    Unobserved stochastic nodes: 59
##    Total graph size: 1520
## 
## Initializing model
 #simulate a sample from the posterior 
#note that we specify both mu and tau variables
crime_sim_1 <- coda.samples(crime_jags_1, variable.names=c("p", "delta"), n.iter=60000)


#store the samples in a data frame:    
crime_data_1 <- data.frame(crime_sim_1[[1]])

Parameter Interpretation

ggplot(crime_data_1, aes(x=delta.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Crime Differential Growth in 2010", x="Crime Differential Growth", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(crime_data_1, aes(x=p.1.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2010 in Zipcode 1", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Confidence Intervals and Posterior Means of Crime Model 1 Parameters

crime_model_1_CI <- matrix(nrow=length(names(crime_data_1)), ncol=2)
for (i in 1:length(names(crime_data_1))) {
      crime_model_1_CI[i, ] <- quantile(crime_data_1[[i]], c(0.025, 0.975))
}

crime_model_1_mean <- matrix(nrow=length(names(crime_data_1)), ncol=1)
for (i in 1:length(names(crime_data_1))) {
      crime_model_1_mean[i, ] <- mean(crime_data_1[[i]])
}

crime_model_1_table <- cbind(crime_model_1_CI, crime_model_1_mean)
rownames(crime_model_1_table) <- names(crime_data_1)
colnames(crime_model_1_table) <- c("2.5%", "97.5%", "Mean")
as.table(crime_model_1_table)
##                    2.5%         97.5%          Mean
## delta.1.  -3.711729e-03 -3.456400e-03 -3.580096e-03
## delta.2.  -2.867138e-03 -2.586290e-03 -2.722396e-03
## delta.3.  -2.969974e-03 -2.586799e-03 -2.769938e-03
## delta.4.  -2.718293e-03 -2.298821e-03 -2.498395e-03
## delta.5.  -4.051228e-03 -3.821873e-03 -3.933484e-03
## delta.6.  -3.588951e-03 -3.197713e-03 -3.382634e-03
## delta.7.  -3.879030e-03 -3.684672e-03 -3.779524e-03
## delta.8.  -4.119015e-03 -3.975270e-03 -4.046257e-03
## delta.9.  -3.750831e-03 -3.637632e-03 -3.693651e-03
## delta.10. -4.088548e-03 -3.894629e-03 -3.989673e-03
## delta.11. -3.756364e-03 -3.595060e-03 -3.674093e-03
## delta.12. -3.611611e-03 -3.482366e-03 -3.545865e-03
## delta.13. -4.320107e-03 -4.105318e-03 -4.209816e-03
## delta.14. -4.097897e-03 -3.946837e-03 -4.021145e-03
## delta.15. -3.980550e-03 -3.808793e-03 -3.892875e-03
## delta.16. -4.056067e-03 -3.886810e-03 -3.970244e-03
## delta.17. -3.837891e-03 -3.730789e-03 -3.783815e-03
## delta.18. -4.166437e-03 -4.030778e-03 -4.097490e-03
## delta.19. -3.668010e-03 -3.564934e-03 -3.615850e-03
## delta.20. -3.587042e-03 -3.497329e-03 -3.541414e-03
## delta.21. -3.400624e-03 -3.292082e-03 -3.346049e-03
## delta.22. -3.974423e-03 -3.823658e-03 -3.897974e-03
## delta.23. -3.749702e-03 -3.652031e-03 -3.700433e-03
## delta.24. -3.351117e-03 -3.254202e-03 -3.302190e-03
## delta.25. -4.264045e-03 -4.100337e-03 -4.180989e-03
## delta.26. -4.122069e-03 -3.945685e-03 -4.032168e-03
## delta.27. -3.631162e-03 -3.536932e-03 -3.583492e-03
## delta.28. -3.966688e-03 -3.860937e-03 -3.913238e-03
## delta.29. -4.278376e-03 -4.086159e-03 -4.179781e-03
## delta.30. -4.747881e-03 -4.358148e-03 -4.543484e-03
## delta.31. -4.150986e-03 -4.012304e-03 -4.080392e-03
## delta.32. -4.613369e-03 -4.136416e-03 -4.360774e-03
## delta.33. -4.320027e-03 -4.140211e-03 -4.228669e-03
## delta.34. -3.480584e-03 -3.371023e-03 -3.425328e-03
## delta.35. -3.604060e-03 -3.491357e-03 -3.546880e-03
## delta.36. -4.438547e-03 -4.213139e-03 -4.321765e-03
## delta.37. -4.011126e-03 -3.890165e-03 -3.949733e-03
## delta.38. -4.046383e-03 -3.898156e-03 -3.970756e-03
## delta.39. -4.120865e-03 -3.969203e-03 -4.043773e-03
## delta.40. -3.825775e-03 -3.611044e-03 -3.715930e-03
## delta.41. -3.908333e-03 -3.764577e-03 -3.835530e-03
## delta.42. -3.463873e-03 -3.365678e-03 -3.414195e-03
## delta.43. -4.320507e-03 -4.100282e-03 -4.207425e-03
## delta.44. -4.819598e-03 -4.399923e-03 -4.598882e-03
## delta.45. -4.029310e-03 -3.903263e-03 -3.965421e-03
## delta.46. -3.567763e-03 -3.456856e-03 -3.511583e-03
## delta.47. -3.645227e-03 -3.541671e-03 -3.592990e-03
## delta.48. -4.309039e-03 -4.080336e-03 -4.191130e-03
## delta.49. -3.715558e-03 -3.564806e-03 -3.638507e-03
## delta.50. -3.663092e-03 -3.465155e-03 -3.561759e-03
## delta.51. -4.471091e-03 -4.157674e-03 -4.307555e-03
## delta.52. -4.593938e-03 -4.240701e-03 -4.410293e-03
## delta.53. -4.154460e-03 -3.995037e-03 -4.072562e-03
## delta.54. -4.328963e-03 -4.090126e-03 -4.206603e-03
## delta.55. -4.244785e-03 -4.029039e-03 -4.134586e-03
## delta.56. -3.888913e-03 -3.554766e-03 -3.714995e-03
## delta.57. -4.794647e-03 -4.448463e-03 -4.614012e-03
## delta.58. -4.593432e-03 -4.245760e-03 -4.413045e-03
## p.1.1.     5.749939e-04  9.602428e-04  7.554488e-04
## p.2.1.     3.131995e-03  5.494800e-03  4.227676e-03
## p.3.1.     2.548631e-03  5.489217e-03  3.877168e-03
## p.4.1.     4.219669e-03  9.750573e-03  6.699385e-03
## p.5.1.     2.906884e-04  4.608558e-04  3.708225e-04
## p.6.1.     7.358186e-04  1.614035e-03  1.136081e-03
## p.7.1.     4.108583e-04  6.071116e-04  5.042648e-04
## p.8.1.     2.536708e-04  3.386202e-04  2.943978e-04
## p.9.1.     5.315555e-04  6.672752e-04  5.972580e-04
## p.10.1.    2.696861e-04  3.981812e-04  3.305846e-04
## p.11.1.    5.256803e-04  7.268444e-04  6.222627e-04
## p.12.1.    7.030785e-04  9.114562e-04  8.040822e-04
## p.13.1.    1.693430e-04  2.607494e-04  2.126423e-04
## p.14.1.    2.646669e-04  3.585291e-04  3.097355e-04
## p.15.1.    3.350469e-04  4.731263e-04  4.011305e-04
## p.16.1.    2.878756e-04  4.044865e-04  3.433503e-04
## p.17.1.    4.462611e-04  5.533940e-04  4.982303e-04
## p.18.1.    2.306134e-04  3.028821e-04  2.655262e-04
## p.19.1.    6.277760e-04  7.721824e-04  6.980831e-04
## p.20.1.    7.386444e-04  8.844757e-04  8.103880e-04
## p.21.1.    1.074039e-03  1.335537e-03  1.200276e-03
## p.22.1.    3.391972e-04  4.592063e-04  3.966882e-04
## p.23.1.    5.327634e-04  6.482515e-04  5.889284e-04
## p.24.1.    1.186280e-03  1.441042e-03  1.310311e-03
## p.25.1.    1.895382e-04  2.633725e-04  2.247577e-04
## p.26.1.    2.521184e-04  3.593598e-04  3.032732e-04
## p.27.1.    6.760041e-04  8.168536e-04  7.448089e-04
## p.28.1.    3.445098e-04  4.260690e-04  3.841264e-04
## p.29.1.    1.841574e-04  2.709842e-04  2.256018e-04
## p.30.1.    7.167902e-05  1.568791e-04  1.102606e-04
## p.31.1.    2.378859e-04  3.143366e-04  2.748372e-04
## p.32.1.    9.392947e-05  2.449535e-04  1.606904e-04
## p.33.1.    1.693701e-04  2.430927e-04  2.043642e-04
## p.34.1.    9.147224e-04  1.139806e-03  1.023656e-03
## p.35.1.    7.138241e-04  8.951459e-04  8.020259e-04
## p.36.1.    1.334728e-04  2.099544e-04  1.699001e-04
## p.37.1.    3.150820e-04  4.017686e-04  3.571354e-04
## p.38.1.    2.935322e-04  3.953701e-04  3.427030e-04
## p.39.1.    2.527294e-04  3.427737e-04  2.959657e-04
## p.40.1.    4.572566e-04  7.038798e-04  5.736094e-04
## p.41.1.    3.873673e-04  5.170781e-04  4.495990e-04
## p.42.1.    9.459398e-04  1.152103e-03  1.046489e-03
## p.43.1.    1.692068e-04  2.634018e-04  2.137256e-04
## p.44.1.    6.205724e-05  1.442461e-04  9.893437e-05
## p.45.1.    3.037772e-04  3.913338e-04  3.461138e-04
## p.46.1.    7.678073e-04  9.593636e-04  8.608960e-04
## p.47.1.    6.571725e-04  8.091164e-04  7.309015e-04
## p.48.1.    1.731517e-04  2.741736e-04  2.209596e-04
## p.49.1.    5.705887e-04  7.723806e-04  6.680907e-04
## p.50.1.    6.340084e-04  9.435087e-04  7.811158e-04
## p.51.1.    1.250224e-04  2.347099e-04  1.759274e-04
## p.52.1.    9.767028e-05  1.986417e-04  1.435782e-04
## p.53.1.    2.362311e-04  3.254346e-04  2.794207e-04
## p.54.1.    1.663558e-04  2.688327e-04  2.143053e-04
## p.55.1.    1.970181e-04  3.039426e-04  2.473618e-04
## p.56.1.    4.027812e-04  7.881140e-04  5.795625e-04
## p.57.1.    6.524867e-05  1.308392e-04  9.527735e-05
## p.58.1.    9.776976e-05  1.966325e-04  1.427305e-04
## p.1.2.     5.728648e-04  9.569328e-04  7.527575e-04
## p.2.2.     3.123056e-03  5.480685e-03  4.216274e-03
## p.3.2.     2.541093e-03  5.475114e-03  3.866558e-03
## p.4.2.     4.208263e-03  9.728402e-03  6.682934e-03
## p.5.2.     2.895134e-04  4.590987e-04  3.693698e-04
## p.6.2.     7.331844e-04  1.608890e-03  1.132271e-03
## p.7.2.     4.092683e-04  6.048801e-04  5.023660e-04
## p.8.2.     2.526283e-04  3.372772e-04  2.932102e-04
## p.9.2.     5.295666e-04  6.648539e-04  5.950583e-04
## p.10.2.    2.685860e-04  3.966341e-04  3.292704e-04
## p.11.2.    5.237104e-04  7.242380e-04  6.199842e-04
## p.12.2.    7.005456e-04  9.082905e-04  8.012402e-04
## p.13.2.    1.686131e-04  2.596814e-04  2.117505e-04
## p.14.2.    2.635849e-04  3.571173e-04  3.084938e-04
## p.15.2.    3.337164e-04  4.713286e-04  3.995741e-04
## p.16.2.    2.867107e-04  4.029180e-04  3.419915e-04
## p.17.2.    4.445524e-04  5.513344e-04  4.963503e-04
## p.18.2.    2.296547e-04  3.016641e-04  2.644414e-04
## p.19.2.    6.254789e-04  7.694366e-04  6.955662e-04
## p.20.2.    7.360015e-04  8.813905e-04  8.075264e-04
## p.21.2.    1.070397e-03  1.331154e-03  1.196273e-03
## p.22.2.    3.378522e-04  4.574546e-04  3.951467e-04
## p.23.2.    5.307705e-04  6.458899e-04  5.867551e-04
## p.24.2.    1.182316e-03  1.436367e-03  1.305999e-03
## p.25.2.    1.887319e-04  2.622950e-04  2.238210e-04
## p.26.2.    2.510815e-04  3.579452e-04  3.020545e-04
## p.27.2.    6.735556e-04  8.139719e-04  7.421475e-04
## p.28.2.    3.431464e-04  4.244278e-04  3.826273e-04
## p.29.2.    1.833713e-04  2.698795e-04  2.246621e-04
## p.30.2.    7.133952e-05  1.561970e-04  1.097630e-04
## p.31.2.    2.369008e-04  3.130783e-04  2.737190e-04
## p.32.2.    9.349718e-05  2.439426e-04  1.599960e-04
## p.33.2.    1.686401e-04  2.420886e-04  2.035028e-04
## p.34.2.    9.115470e-04  1.135974e-03  1.020161e-03
## p.35.2.    7.112579e-04  8.920289e-04  7.991898e-04
## p.36.2.    1.328817e-04  2.090719e-04  1.691687e-04
## p.37.2.    3.138211e-04  4.002093e-04  3.557287e-04
## p.38.2.    2.923472e-04  3.938325e-04  3.413463e-04
## p.39.2.    2.516904e-04  3.414163e-04  2.947725e-04
## p.40.2.    4.555114e-04  7.013444e-04  5.714865e-04
## p.41.2.    3.858569e-04  5.151362e-04  4.478798e-04
## p.42.2.    9.426719e-04  1.148236e-03  1.042927e-03
## p.43.2.    1.684774e-04  2.623243e-04  2.128298e-04
## p.44.2.    6.175889e-05  1.436129e-04  9.848270e-05
## p.45.2.    3.025560e-04  3.898099e-04  3.447452e-04
## p.46.2.    7.650749e-04  9.560561e-04  8.578822e-04
## p.47.2.    6.547829e-04  8.062581e-04  7.282830e-04
## p.48.2.    1.724074e-04  2.730575e-04  2.200372e-04
## p.49.2.    5.684738e-04  7.696343e-04  6.656679e-04
## p.50.2.    6.316917e-04  9.402480e-04  7.783447e-04
## p.51.2.    1.244647e-04  2.337363e-04  1.751736e-04
## p.52.2.    9.722266e-05  1.978012e-04  1.429488e-04
## p.53.2.    2.352520e-04  3.241375e-04  2.782863e-04
## p.54.2.    1.656373e-04  2.677357e-04  2.134074e-04
## p.55.2.    1.961838e-04  3.027209e-04  2.463430e-04
## p.56.2.    4.012185e-04  7.853196e-04  5.774230e-04
## p.57.2.    6.493660e-05  1.302585e-04  9.484026e-05
## p.58.2.    9.732173e-05  1.957996e-04  1.421043e-04
## p.1.3.     5.707437e-04  9.536341e-04  7.500757e-04
## p.2.3.     3.114142e-03  5.466606e-03  4.204903e-03
## p.3.3.     2.533576e-03  5.461046e-03  3.855977e-03
## p.4.3.     4.196887e-03  9.706281e-03  6.666523e-03
## p.5.3.     2.883433e-04  4.573482e-04  3.679228e-04
## p.6.3.     7.305597e-04  1.603762e-03  1.128474e-03
## p.7.3.     4.076845e-04  6.026567e-04  5.004743e-04
## p.8.3.     2.515901e-04  3.359396e-04  2.920273e-04
## p.9.3.     5.275850e-04  6.624414e-04  5.928667e-04
## p.10.3.    2.674904e-04  3.950929e-04  3.279613e-04
## p.11.3.    5.217478e-04  7.216408e-04  6.177140e-04
## p.12.3.    6.980218e-04  9.051359e-04  7.984081e-04
## p.13.3.    1.678864e-04  2.586178e-04  2.108624e-04
## p.14.3.    2.625072e-04  3.557111e-04  3.072571e-04
## p.15.3.    3.323911e-04  4.695376e-04  3.980238e-04
## p.16.3.    2.855504e-04  4.013556e-04  3.406382e-04
## p.17.3.    4.428503e-04  5.492824e-04  4.944775e-04
## p.18.3.    2.287001e-04  3.004510e-04  2.633610e-04
## p.19.3.    6.231903e-04  7.667006e-04  6.930584e-04
## p.20.3.    7.333681e-04  8.783161e-04  8.046748e-04
## p.21.3.    1.066767e-03  1.326784e-03  1.192284e-03
## p.22.3.    3.365125e-04  4.557096e-04  3.936112e-04
## p.23.3.    5.287850e-04  6.435369e-04  5.845899e-04
## p.24.3.    1.178365e-03  1.431707e-03  1.301700e-03
## p.25.3.    1.879290e-04  2.612220e-04  2.228881e-04
## p.26.3.    2.500490e-04  3.565361e-04  3.008406e-04
## p.27.3.    6.711158e-04  8.111003e-04  7.394956e-04
## p.28.3.    3.417884e-04  4.227930e-04  3.811340e-04
## p.29.3.    1.825886e-04  2.687793e-04  2.237263e-04
## p.30.3.    7.100164e-05  1.555179e-04  1.092676e-04
## p.31.3.    2.359196e-04  3.118250e-04  2.726054e-04
## p.32.3.    9.306688e-05  2.429359e-04  1.593045e-04
## p.33.3.    1.679133e-04  2.410886e-04  2.026451e-04
## p.34.3.    9.083827e-04  1.132156e-03  1.016678e-03
## p.35.3.    7.087009e-04  8.889227e-04  7.963638e-04
## p.36.3.    1.322933e-04  2.081931e-04  1.684404e-04
## p.37.3.    3.125652e-04  3.986560e-04  3.543277e-04
## p.38.3.    2.911670e-04  3.923009e-04  3.399951e-04
## p.39.3.    2.506556e-04  3.400643e-04  2.935842e-04
## p.40.3.    4.537728e-04  6.988182e-04  5.693715e-04
## p.41.3.    3.843524e-04  5.132016e-04  4.461672e-04
## p.42.3.    9.394154e-04  1.144383e-03  1.039378e-03
## p.43.3.    1.677512e-04  2.612512e-04  2.119377e-04
## p.44.3.    6.146197e-05  1.429825e-04  9.803308e-05
## p.45.3.    3.013397e-04  3.882920e-04  3.433821e-04
## p.46.3.    7.623523e-04  9.527600e-04  8.548789e-04
## p.47.3.    6.524019e-04  8.034100e-04  7.256739e-04
## p.48.3.    1.716662e-04  2.719459e-04  2.191186e-04
## p.49.3.    5.663667e-04  7.668976e-04  6.632538e-04
## p.50.3.    6.293834e-04  9.369986e-04  7.755835e-04
## p.51.3.    1.239095e-04  2.327668e-04  1.744230e-04
## p.52.3.    9.677709e-05  1.969644e-04  1.423221e-04
## p.53.3.    2.342769e-04  3.228456e-04  2.771565e-04
## p.54.3.    1.649220e-04  2.666432e-04  2.125134e-04
## p.55.3.    1.953529e-04  3.015040e-04  2.453283e-04
## p.56.3.    3.996618e-04  7.825351e-04  5.752914e-04
## p.57.3.    6.462602e-05  1.296804e-04  9.440517e-05
## p.58.3.    9.687576e-05  1.949702e-04  1.414809e-04
## p.1.4.     5.686303e-04  9.503467e-04  7.474036e-04
## p.2.4.     3.105254e-03  5.452563e-03  4.193562e-03
## p.3.4.     2.526081e-03  5.447015e-03  3.845425e-03
## p.4.4.     4.185542e-03  9.684209e-03  6.650153e-03
## p.5.4.     2.871778e-04  4.556044e-04  3.664815e-04
## p.6.4.     7.279443e-04  1.598650e-03  1.124690e-03
## p.7.4.     4.061068e-04  6.004415e-04  4.985897e-04
## p.8.4.     2.505562e-04  3.346072e-04  2.908492e-04
## p.9.4.     5.256109e-04  6.600377e-04  5.906832e-04
## p.10.4.    2.663993e-04  3.935578e-04  3.266575e-04
## p.11.4.    5.197927e-04  7.190530e-04  6.154521e-04
## p.12.4.    6.955071e-04  9.019922e-04  7.955861e-04
## p.13.4.    1.671628e-04  2.575585e-04  2.099780e-04
## p.14.4.    2.614340e-04  3.543105e-04  3.060253e-04
## p.15.4.    3.310710e-04  4.677535e-04  3.964795e-04
## p.16.4.    2.843949e-04  3.997992e-04  3.392902e-04
## p.17.4.    4.411547e-04  5.472381e-04  4.926116e-04
## p.18.4.    2.277494e-04  2.992427e-04  2.622850e-04
## p.19.4.    6.209100e-04  7.639744e-04  6.905596e-04
## p.20.4.    7.307441e-04  8.752524e-04  8.018333e-04
## p.21.4.    1.063149e-03  1.322429e-03  1.188308e-03
## p.22.4.    3.351782e-04  4.539713e-04  3.920817e-04
## p.23.4.    5.268070e-04  6.411925e-04  5.824326e-04
## p.24.4.    1.174427e-03  1.427062e-03  1.297416e-03
## p.25.4.    1.871295e-04  2.601534e-04  2.219592e-04
## p.26.4.    2.490206e-04  3.551326e-04  2.996316e-04
## p.27.4.    6.686850e-04  8.082389e-04  7.368532e-04
## p.28.4.    3.404358e-04  4.211644e-04  3.796466e-04
## p.29.4.    1.818092e-04  2.676835e-04  2.227944e-04
## p.30.4.    7.066536e-05  1.548417e-04  1.087745e-04
## p.31.4.    2.349426e-04  3.105768e-04  2.714963e-04
## p.32.4.    9.263855e-05  2.419334e-04  1.586160e-04
## p.33.4.    1.671896e-04  2.400928e-04  2.017910e-04
## p.34.4.    9.052294e-04  1.128350e-03  1.013206e-03
## p.35.4.    7.061531e-04  8.858273e-04  7.935478e-04
## p.36.4.    1.317075e-04  2.073180e-04  1.677152e-04
## p.37.4.    3.113144e-04  3.971088e-04  3.529321e-04
## p.38.4.    2.899915e-04  3.907752e-04  3.386492e-04
## p.39.4.    2.496250e-04  3.387177e-04  2.924006e-04
## p.40.4.    4.520409e-04  6.963010e-04  5.672643e-04
## p.41.4.    3.828537e-04  5.112742e-04  4.444612e-04
## p.42.4.    9.361700e-04  1.140542e-03  1.035840e-03
## p.43.4.    1.670281e-04  2.601825e-04  2.110494e-04
## p.44.4.    6.116648e-05  1.423549e-04  9.758552e-05
## p.45.4.    3.001283e-04  3.867799e-04  3.420243e-04
## p.46.4.    7.596393e-04  9.494752e-04  8.518861e-04
## p.47.4.    6.500296e-04  8.005719e-04  7.230742e-04
## p.48.4.    1.709282e-04  2.708388e-04  2.182039e-04
## p.49.4.    5.642674e-04  7.641708e-04  6.608485e-04
## p.50.4.    6.270836e-04  9.337604e-04  7.728321e-04
## p.51.4.    1.233568e-04  2.318012e-04  1.736756e-04
## p.52.4.    9.633356e-05  1.961310e-04  1.416982e-04
## p.53.4.    2.333058e-04  3.215588e-04  2.760313e-04
## p.54.4.    1.642097e-04  2.655551e-04  2.116230e-04
## p.55.4.    1.945256e-04  3.002921e-04  2.443178e-04
## p.56.4.    3.981112e-04  7.797605e-04  5.731676e-04
## p.57.4.    6.431692e-05  1.291049e-04  9.397207e-05
## p.58.4.    9.643183e-05  1.941443e-04  1.408602e-04
## p.1.5.     5.665248e-04  9.470707e-04  7.447409e-04
## p.2.5.     3.096391e-03  5.438556e-03  4.182252e-03
## p.3.5.     2.518609e-03  5.433019e-03  3.834902e-03
## p.4.5.     4.174227e-03  9.662187e-03  6.633822e-03
## p.5.5.     2.860171e-04  4.538672e-04  3.650458e-04
## p.6.5.     7.253384e-04  1.593554e-03  1.120919e-03
## p.7.5.     4.045352e-04  5.982345e-04  4.967122e-04
## p.8.5.     2.495265e-04  3.332802e-04  2.896758e-04
## p.9.5.     5.236441e-04  6.576426e-04  5.885077e-04
## p.10.5.    2.653126e-04  3.920286e-04  3.253589e-04
## p.11.5.    5.178448e-04  7.164745e-04  6.131985e-04
## p.12.5.    6.930015e-04  8.988594e-04  7.927740e-04
## p.13.5.    1.664423e-04  2.565036e-04  2.090973e-04
## p.14.5.    2.603651e-04  3.529153e-04  3.047985e-04
## p.15.5.    3.297563e-04  4.659761e-04  3.949412e-04
## p.16.5.    2.832440e-04  3.982489e-04  3.379475e-04
## p.17.5.    4.394656e-04  5.452014e-04  4.907529e-04
## p.18.5.    2.268027e-04  2.980393e-04  2.612134e-04
## p.19.5.    6.186381e-04  7.612578e-04  6.880698e-04
## p.20.5.    7.281295e-04  8.721994e-04  7.990018e-04
## p.21.5.    1.059544e-03  1.318089e-03  1.184345e-03
## p.22.5.    3.338492e-04  4.522395e-04  3.905581e-04
## p.23.5.    5.248364e-04  6.388566e-04  5.802833e-04
## p.24.5.    1.170503e-03  1.422432e-03  1.293146e-03
## p.25.5.    1.863334e-04  2.590891e-04  2.210341e-04
## p.26.5.    2.479965e-04  3.537346e-04  2.984274e-04
## p.27.5.    6.662629e-04  8.053876e-04  7.342202e-04
## p.28.5.    3.390885e-04  4.195422e-04  3.781650e-04
## p.29.5.    1.810332e-04  2.665922e-04  2.218664e-04
## p.30.5.    7.033066e-05  1.541684e-04  1.082836e-04
## p.31.5.    2.339696e-04  3.093335e-04  2.703917e-04
## p.32.5.    9.221220e-05  2.409349e-04  1.579305e-04
## p.33.5.    1.664690e-04  2.391010e-04  2.009406e-04
## p.34.5.    9.020869e-04  1.124557e-03  1.009747e-03
## p.35.5.    7.036145e-04  8.827427e-04  7.907417e-04
## p.36.5.    1.311243e-04  2.064465e-04  1.669932e-04
## p.37.5.    3.100685e-04  3.955676e-04  3.515420e-04
## p.38.5.    2.888208e-04  3.892555e-04  3.373086e-04
## p.39.5.    2.485987e-04  3.373763e-04  2.912218e-04
## p.40.5.    4.503156e-04  6.937929e-04  5.651649e-04
## p.41.5.    3.813609e-04  5.093541e-04  4.427617e-04
## p.42.5.    9.329359e-04  1.136714e-03  1.032314e-03
## p.43.5.    1.663081e-04  2.591181e-04  2.101648e-04
## p.44.5.    6.087241e-05  1.417300e-04  9.714000e-05
## p.45.5.    2.989218e-04  3.852737e-04  3.406719e-04
## p.46.5.    7.569359e-04  9.462018e-04  8.489037e-04
## p.47.5.    6.476660e-04  7.977438e-04  7.204837e-04
## p.48.5.    1.701933e-04  2.697363e-04  2.172930e-04
## p.49.5.    5.621759e-04  7.614536e-04  6.584519e-04
## p.50.5.    6.247921e-04  9.305334e-04  7.700904e-04
## p.51.5.    1.228066e-04  2.308397e-04  1.729315e-04
## p.52.5.    9.589207e-05  1.953012e-04  1.410770e-04
## p.53.5.    2.323388e-04  3.202771e-04  2.749106e-04
## p.54.5.    1.635005e-04  2.644714e-04  2.107364e-04
## p.55.5.    1.937018e-04  2.990850e-04  2.433114e-04
## p.56.5.    3.965666e-04  7.769957e-04  5.710517e-04
## p.57.5.    6.400930e-05  1.285319e-04  9.354097e-05
## p.58.5.    9.598994e-05  1.933219e-04  1.402422e-04
## p.1.6.     5.644271e-04  9.438060e-04  7.420878e-04
## p.2.6.     3.087553e-03  5.424585e-03  4.170972e-03
## p.3.6.     2.511159e-03  5.419059e-03  3.824407e-03
## p.4.6.     4.162943e-03  9.640215e-03  6.617532e-03
## p.5.6.     2.848610e-04  4.521367e-04  3.636157e-04
## p.6.6.     7.227417e-04  1.588474e-03  1.117160e-03
## p.7.6.     4.029696e-04  5.960356e-04  4.948418e-04
## p.8.6.     2.485011e-04  3.319584e-04  2.885072e-04
## p.9.6.     5.216847e-04  6.552563e-04  5.863403e-04
## p.10.6.    2.642304e-04  3.905054e-04  3.240654e-04
## p.11.6.    5.159042e-04  7.139052e-04  6.109531e-04
## p.12.6.    6.905049e-04  8.957375e-04  7.899718e-04
## p.13.6.    1.657249e-04  2.554530e-04  2.082204e-04
## p.14.6.    2.593006e-04  3.515256e-04  3.035766e-04
## p.15.6.    3.284467e-04  4.642055e-04  3.934088e-04
## p.16.6.    2.820978e-04  3.967046e-04  3.366102e-04
## p.17.6.    4.377829e-04  5.431723e-04  4.889011e-04
## p.18.6.    2.258599e-04  2.968408e-04  2.601462e-04
## p.19.6.    6.163745e-04  7.585508e-04  6.855890e-04
## p.20.6.    7.255242e-04  8.691570e-04  7.961803e-04
## p.21.6.    1.055950e-03  1.313762e-03  1.180395e-03
## p.22.6.    3.325254e-04  4.505144e-04  3.890404e-04
## p.23.6.    5.228731e-04  6.365292e-04  5.781420e-04
## p.24.6.    1.166592e-03  1.417818e-03  1.288890e-03
## p.25.6.    1.855407e-04  2.580292e-04  2.201129e-04
## p.26.6.    2.469766e-04  3.523422e-04  2.972281e-04
## p.27.6.    6.638495e-04  8.025463e-04  7.315966e-04
## p.28.6.    3.377466e-04  4.179261e-04  3.766891e-04
## p.29.6.    1.802604e-04  2.655054e-04  2.209423e-04
## p.30.6.    6.999756e-05  1.534981e-04  1.077949e-04
## p.31.6.    2.330006e-04  3.080953e-04  2.692917e-04
## p.32.6.    9.178781e-05  2.399406e-04  1.572480e-04
## p.33.6.    1.657515e-04  2.381134e-04  2.000937e-04
## p.34.6.    8.989554e-04  1.120777e-03  1.006299e-03
## p.35.6.    7.010849e-04  8.796688e-04  7.879455e-04
## p.36.6.    1.305436e-04  2.055787e-04  1.662743e-04
## p.37.6.    3.088277e-04  3.940324e-04  3.501574e-04
## p.38.6.    2.876548e-04  3.877416e-04  3.359733e-04
## p.39.6.    2.475766e-04  3.360403e-04  2.900478e-04
## p.40.6.    4.485968e-04  6.912939e-04  5.630733e-04
## p.41.6.    3.798738e-04  5.074411e-04  4.410686e-04
## p.42.6.    9.297129e-04  1.132899e-03  1.028801e-03
## p.43.6.    1.655913e-04  2.580581e-04  2.092839e-04
## p.44.6.    6.057975e-05  1.411078e-04  9.669652e-05
## p.45.6.    2.977201e-04  3.837734e-04  3.393248e-04
## p.46.6.    7.542422e-04  9.429396e-04  8.459319e-04
## p.47.6.    6.453109e-04  7.949257e-04  7.179025e-04
## p.48.6.    1.694617e-04  2.686382e-04  2.163858e-04
## p.49.6.    5.600921e-04  7.587460e-04  6.560640e-04
## p.50.6.    6.225091e-04  9.273175e-04  7.673584e-04
## p.51.6.    1.222588e-04  2.298822e-04  1.721905e-04
## p.52.6.    9.545260e-05  1.944749e-04  1.404585e-04
## p.53.6.    2.313758e-04  3.190005e-04  2.737945e-04
## p.54.6.    1.627943e-04  2.633922e-04  2.098535e-04
## p.55.6.    1.928815e-04  2.978827e-04  2.423092e-04
## p.56.6.    3.950280e-04  7.742407e-04  5.689436e-04
## p.57.6.    6.370315e-05  1.279615e-04  9.311184e-05
## p.58.6.    9.555007e-05  1.925030e-04  1.396270e-04
## p.1.7.     5.623372e-04  9.405525e-04  7.394441e-04
## p.2.7.     3.078741e-03  5.410649e-03  4.159722e-03
## p.3.7.     2.503730e-03  5.405135e-03  3.813941e-03
## p.4.7.     4.151689e-03  9.618292e-03  6.601281e-03
## p.5.7.     2.837097e-04  4.504128e-04  3.621913e-04
## p.6.7.     7.201543e-04  1.583411e-03  1.113414e-03
## p.7.7.     4.014102e-04  5.938447e-04  4.929784e-04
## p.8.7.     2.474799e-04  3.306418e-04  2.873433e-04
## p.9.7.     5.197326e-04  6.528786e-04  5.841808e-04
## p.10.7.    2.631525e-04  3.889881e-04  3.227771e-04
## p.11.7.    5.139709e-04  7.113451e-04  6.087160e-04
## p.12.7.    6.880173e-04  8.926264e-04  7.871796e-04
## p.13.7.    1.650107e-04  2.544067e-04  2.073471e-04
## p.14.7.    2.582405e-04  3.501414e-04  3.023597e-04
## p.15.7.    3.271423e-04  4.624417e-04  3.918824e-04
## p.16.7.    2.809563e-04  3.951663e-04  3.352781e-04
## p.17.7.    4.361067e-04  5.411507e-04  4.870563e-04
## p.18.7.    2.249211e-04  2.956470e-04  2.590833e-04
## p.19.7.    6.141192e-04  7.558535e-04  6.831171e-04
## p.20.7.    7.229283e-04  8.661252e-04  7.933687e-04
## p.21.7.    1.052369e-03  1.309450e-03  1.176459e-03
## p.22.7.    3.312068e-04  4.487958e-04  3.875286e-04
## p.23.7.    5.209172e-04  6.342103e-04  5.760085e-04
## p.24.7.    1.162693e-03  1.413218e-03  1.284648e-03
## p.25.7.    1.847514e-04  2.569737e-04  2.191955e-04
## p.26.7.    2.459609e-04  3.509551e-04  2.960336e-04
## p.27.7.    6.614450e-04  7.997150e-04  7.289824e-04
## p.28.7.    3.364099e-04  4.163163e-04  3.752190e-04
## p.29.7.    1.794910e-04  2.644230e-04  2.200220e-04
## p.30.7.    6.966603e-05  1.528307e-04  1.073084e-04
## p.31.7.    2.320357e-04  3.068619e-04  2.681961e-04
## p.32.7.    9.136537e-05  2.389504e-04  1.565684e-04
## p.33.7.    1.650371e-04  2.371298e-04  1.992503e-04
## p.34.7.    8.958348e-04  1.117009e-03  1.002863e-03
## p.35.7.    6.985645e-04  8.766056e-04  7.851592e-04
## p.36.7.    1.299656e-04  2.047146e-04  1.655584e-04
## p.37.7.    3.075918e-04  3.925031e-04  3.487783e-04
## p.38.7.    2.864935e-04  3.862337e-04  3.346433e-04
## p.39.7.    2.465588e-04  3.347096e-04  2.888784e-04
## p.40.7.    4.468846e-04  6.888038e-04  5.609893e-04
## p.41.7.    3.783926e-04  5.055354e-04  4.393821e-04
## p.42.7.    9.265010e-04  1.129097e-03  1.025299e-03
## p.43.7.    1.648775e-04  2.570024e-04  2.084067e-04
## p.44.7.    6.028850e-05  1.404884e-04  9.625506e-05
## p.45.7.    2.965233e-04  3.822789e-04  3.379831e-04
## p.46.7.    7.515581e-04  9.396887e-04  8.429704e-04
## p.47.7.    6.429644e-04  7.921175e-04  7.153306e-04
## p.48.7.    1.687332e-04  2.675446e-04  2.154825e-04
## p.49.7.    5.580161e-04  7.560481e-04  6.536847e-04
## p.50.7.    6.202344e-04  9.241127e-04  7.646362e-04
## p.51.7.    1.217134e-04  2.289286e-04  1.714527e-04
## p.52.7.    9.501515e-05  1.936521e-04  1.398428e-04
## p.53.7.    2.304168e-04  3.177291e-04  2.726830e-04
## p.54.7.    1.620912e-04  2.623174e-04  2.089743e-04
## p.55.7.    1.920646e-04  2.966853e-04  2.413112e-04
## p.56.7.    3.934953e-04  7.714955e-04  5.668433e-04
## p.57.7.    6.339847e-05  1.273936e-04  9.268468e-05
## p.58.7.    9.511221e-05  1.916876e-04  1.390144e-04
## p.1.8.     5.602550e-04  9.373103e-04  7.368098e-04
## p.2.8.     3.069953e-03  5.396749e-03  4.148503e-03
## p.3.8.     2.496324e-03  5.391246e-03  3.803504e-03
## p.4.8.     4.140466e-03  9.596419e-03  6.585070e-03
## p.5.8.     2.825629e-04  4.486954e-04  3.607724e-04
## p.6.8.     7.175762e-04  1.578364e-03  1.109681e-03
## p.7.8.     3.998567e-04  5.916619e-04  4.911220e-04
## p.8.8.     2.464629e-04  3.293304e-04  2.861841e-04
## p.9.8.     5.177879e-04  6.505095e-04  5.820292e-04
## p.10.8.    2.620791e-04  3.874766e-04  3.214939e-04
## p.11.8.    5.120449e-04  7.087941e-04  6.064870e-04
## p.12.8.    6.855386e-04  8.895261e-04  7.843972e-04
## p.13.8.    1.642994e-04  2.533647e-04  2.064774e-04
## p.14.8.    2.571847e-04  3.487627e-04  3.011475e-04
## p.15.8.    3.258431e-04  4.606845e-04  3.903620e-04
## p.16.8.    2.798193e-04  3.936340e-04  3.339513e-04
## p.17.8.    4.344369e-04  5.391366e-04  4.852185e-04
## p.18.8.    2.239861e-04  2.944581e-04  2.580248e-04
## p.19.8.    6.118721e-04  7.531657e-04  6.806541e-04
## p.20.8.    7.203416e-04  8.631040e-04  7.905671e-04
## p.21.8.    1.048800e-03  1.305152e-03  1.172535e-03
## p.22.8.    3.298935e-04  4.470838e-04  3.860227e-04
## p.23.8.    5.189685e-04  6.318999e-04  5.738829e-04
## p.24.8.    1.158808e-03  1.408633e-03  1.280420e-03
## p.25.8.    1.839655e-04  2.559224e-04  2.182819e-04
## p.26.8.    2.449494e-04  3.495736e-04  2.948439e-04
## p.27.8.    6.590491e-04  7.968937e-04  7.263775e-04
## p.28.8.    3.350786e-04  4.147127e-04  3.737547e-04
## p.29.8.    1.787248e-04  2.633450e-04  2.191055e-04
## p.30.8.    6.933607e-05  1.521662e-04  1.068241e-04
## p.31.8.    2.310747e-04  3.056336e-04  2.671049e-04
## p.32.8.    9.094488e-05  2.379643e-04  1.558917e-04
## p.33.8.    1.643258e-04  2.361503e-04  1.984105e-04
## p.34.8.    8.927250e-04  1.113254e-03  9.994387e-04
## p.35.8.    6.960531e-04  8.735531e-04  7.823827e-04
## p.36.8.    1.293901e-04  2.038541e-04  1.648457e-04
## p.37.8.    3.063609e-04  3.909798e-04  3.474046e-04
## p.38.8.    2.853370e-04  3.847316e-04  3.333185e-04
## p.39.8.    2.455451e-04  3.333841e-04  2.877138e-04
## p.40.8.    4.451790e-04  6.863227e-04  5.589132e-04
## p.41.8.    3.769172e-04  5.036368e-04  4.377020e-04
## p.42.8.    9.233003e-04  1.125307e-03  1.021809e-03
## p.43.8.    1.641668e-04  2.559511e-04  2.075332e-04
## p.44.8.    5.999865e-05  1.398717e-04  9.581562e-05
## p.45.8.    2.953313e-04  3.807903e-04  3.366466e-04
## p.46.8.    7.488835e-04  9.364490e-04  8.400192e-04
## p.47.8.    6.406264e-04  7.893193e-04  7.127679e-04
## p.48.8.    1.680078e-04  2.664554e-04  2.145829e-04
## p.49.8.    5.559478e-04  7.533597e-04  6.513141e-04
## p.50.8.    6.179679e-04  9.209190e-04  7.619236e-04
## p.51.8.    1.211705e-04  2.279790e-04  1.707181e-04
## p.52.8.    9.457969e-05  1.928328e-04  1.392297e-04
## p.53.8.    2.294617e-04  3.164627e-04  2.715759e-04
## p.54.8.    1.613912e-04  2.612469e-04  2.080988e-04
## p.55.8.    1.912512e-04  2.954927e-04  2.403172e-04
## p.56.8.    3.919686e-04  7.687600e-04  5.647507e-04
## p.57.8.    6.309524e-05  1.268282e-04  9.225948e-05
## p.58.8.    9.467636e-05  1.908756e-04  1.384045e-04

The 95% confidence interval for \(\beta\) is (-0.004096, 0.0002897). Since \(\beta\) is not considered signifant in the 5% level, we can proceed with Model 2 instead, in which the overall time trend is removed from the model.

Crime Model 1 Visualizations

#Alter the data frame for visualizations
zip_data <- data.frame(zipcode=c(1:58), zip=rownames(Crime_count))

crime_viz_1 <- crime_data_1 %>%
  colMeans() %>%
  as.data.frame() %>%
  rownames_to_column(var="zipcode") %>%
  setNames(., c("zipcode", "delta")) %>%
  filter(grepl("delta", zipcode)) %>%
  mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
  left_join(zip_data, by="zipcode") %>%
  mutate(id=as.character(zip))

Map of Crime Differential Growth by Zipcode (Using Model 1)

zip <- readOGR(dsn="Boundaries - ZIP Codes", layer="geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a")
## OGR data source with driver: ESRI Shapefile 
## Source: "Boundaries - ZIP Codes", layer: "geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a"
## with 61 features
## It has 4 fields
zip_df <- fortify(zip, region = "zip")
crime_model_1_delta <- left_join(zip_df, crime_viz_1, by="id")
head(crime_model_1_delta)
##        long      lat order  hole piece    id   group zipcode        delta
## 1 -87.62271 41.88884     1 FALSE     1 60601 60601.1       1 -0.003580096
## 2 -87.62232 41.88879     2 FALSE     1 60601 60601.1       1 -0.003580096
## 3 -87.62120 41.88866     3 FALSE     1 60601 60601.1       1 -0.003580096
## 4 -87.62064 41.88860     4 FALSE     1 60601 60601.1       1 -0.003580096
## 5 -87.62058 41.88859     5 FALSE     1 60601 60601.1       1 -0.003580096
## 6 -87.62055 41.88859     6 FALSE     1 60601 60601.1       1 -0.003580096
##     zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
model_1_map <- ggplot() 
model_1_map <- model_1_map +  geom_polygon(data = crime_model_1_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2) 
model_1_map <- model_1_map + coord_map() 
model_1_map <- model_1_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_1_map <- model_1_map + cleanup
print(model_1_map)


Model 3 – Use demographics information

Notation

  • i: zipcodes in Chicago (58 total in our data)

  • j: years (from 2010 to 2017)

  • y[i, j]: observed count of violent crimes in Chicago by year and zipcode

  • p[i, j]: probability of being a victim in location i and time j

  • n[i, j]: the population in zipcode i in time j.

  • time[j]: the time in years

  • \(\delta\)[i]: the spatial temporal interaction term

  • \(\sigma\)[i]: the income trend

  • income[i]: income for each zipcode in Chicago

Specify the Model

This model resembles Model 2 but now our prediction of the crime rate accounts for the different income levels in various Chicago zipcodes.

Data

\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \delta_{i} * time_{j} + \sigma_{i} * income_{i}\]

Priors

\[\delta_{i} \sim N(0, prec.delta)\] \[\sigma_{i} \sim N(0, prec.sigma)\]

Hyperprior

\[prec.u \sim Gamma(0.5, 0.0005)\] \[prec.delta \sim Gamma(0.5, 0.0005)\] \[prec.sigma \sim Gamma(0.5, 0.0005)\]

#Specify the model
crime_model_3 <- "model{
    #Data
    for(i in 1:58) {
      for(j in 1:8){
        y[i,j] ~ dbin(p[i,j], n[i])
        logit(p[i,j]) = delta[i] * time[j] + sigma[i] * income[i]
      }
    }
    
    #Priors
    for(i in 1:58){
      delta[i] ~ dnorm(0, prec.delta)
      sigma[i] ~ dnorm(0, prec.sigma)
    }

    #Hyperpriors
    prec.delta ~ dgamma(0.5, 0.0005)
    prec.sigma ~ dgamma(0.5, 0.0005)
    
}"

#set up an algorithm to simulate the posterior by 
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_3 <- jags.model(textConnection(crime_model_3), data=list(y=Crime_count, n=Population.vector, time=time, income=Model_Data$per_capita_income),
                    inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 464
##    Unobserved stochastic nodes: 118
##    Total graph size: 2564
## 
## Initializing model
#simulate a sample from the posterior 
#note that we specify both mu and tau variables
crime_sim_3 <- coda.samples(crime_jags_3, variable.names=c("p", "delta", "sigma"), n.iter=40000)


#store the samples in a data frame:    
crime_data_3 <- data.frame(crime_sim_3[[1]])

Crime Model 3 Visualizations

crime_viz_3_delta <- crime_data_3 %>%
  colMeans() %>%
  as.data.frame() %>%
  rownames_to_column(var="zipcode") %>%
  setNames(., c("zipcode", "delta")) %>%
  filter(grepl("delta", zipcode)) %>%
  mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
  left_join(zip_data, by="zipcode") %>%
  mutate(id=as.character(zip))
crime_model_3_delta <- left_join(zip_df, crime_viz_3_delta, by="id")
head(crime_model_3_delta)
##        long      lat order  hole piece    id   group zipcode      delta
## 1 -87.62271 41.88884     1 FALSE     1 60601 60601.1       1 -0.1297566
## 2 -87.62232 41.88879     2 FALSE     1 60601 60601.1       1 -0.1297566
## 3 -87.62120 41.88866     3 FALSE     1 60601 60601.1       1 -0.1297566
## 4 -87.62064 41.88860     4 FALSE     1 60601 60601.1       1 -0.1297566
## 5 -87.62058 41.88859     5 FALSE     1 60601 60601.1       1 -0.1297566
## 6 -87.62055 41.88859     6 FALSE     1 60601 60601.1       1 -0.1297566
##     zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
crime_viz_3_sigma <- crime_data_3 %>%
  colMeans() %>%
  as.data.frame() %>%
  rownames_to_column(var="zipcode") %>%
  setNames(., c("zipcode", "sigma")) %>%
  filter(grepl("sigma", zipcode)) %>%
  mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
  left_join(zip_data, by="zipcode") %>%
  mutate(id=as.character(zip))
crime_model_3_sigma <- left_join(zip_df, crime_viz_3_sigma, by="id")
head(crime_model_3_sigma)
##        long      lat order  hole piece    id   group zipcode       sigma
## 1 -87.62271 41.88884     1 FALSE     1 60601 60601.1       1 0.002849706
## 2 -87.62232 41.88879     2 FALSE     1 60601 60601.1       1 0.002849706
## 3 -87.62120 41.88866     3 FALSE     1 60601 60601.1       1 0.002849706
## 4 -87.62064 41.88860     4 FALSE     1 60601 60601.1       1 0.002849706
## 5 -87.62058 41.88859     5 FALSE     1 60601 60601.1       1 0.002849706
## 6 -87.62055 41.88859     6 FALSE     1 60601 60601.1       1 0.002849706
##     zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601

Map of Crime Differential Growth by Zipcode (Using Model 3)

model_3_map <- ggplot() 
model_3_map <- model_3_map +  geom_polygon(data = crime_model_3_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2) 
model_3_map <- model_3_map + coord_map() 
model_3_map <- model_3_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_3_map <- model_3_map + cleanup
print(model_3_map)

Income Differential Growth by Zipcode (Using Model 3)

model_3_map <- ggplot() 
model_3_map <- model_3_map +  geom_polygon(data = crime_model_3_sigma, aes(x=long, y=lat, group=group, fill=sigma), color = "black", size=0.2) 
model_3_map <- model_3_map + coord_map() 
model_3_map <- model_3_map + labs(x="", y="", title=paste("Income Differential Growth by Zipcode"), fill="total")
model_3_map <- model_3_map + cleanup
print(model_3_map)

In locations where the income differential growth increases, the crime differential growth seems to decrease.

Running mean plot of alpha

The running mean of \(\alpha\) appears to stabilize after 40,000 iterations.

running_mean_plot(crime_data_3$delta.1.)

running_mean_plot(crime_data_3$delta.3.)

running_mean_plot(crime_data_3$sigma.1.)

running_mean_plot(crime_data_3$sigma.3.)


Model 4 – Use demographics information with race

Notation

  • i: zipcodes in Chicago (58 total in our data)

  • j: years (from 2010 to 2017)

  • y[i, j]: observed count of violent crimes in Chicago by year and zipcode

  • p[i, j]: probability of being a victim in location i and time j

  • n[i, j]: the population in zipcode i in time j.

  • \(\alpha\): the mean crime rate

  • u[i]: the unstructured spatial random effect for each zipcode in Chicago.

  • s[i]: the structured spatial random effect for each zipcode in Chicago.

  • time[j]: the time in years

  • \(\delta\)[i]: the spatial temporal interaction term

  • \(\sigma\)[i]: the income trend

  • income[i]: income for each zipcode in Chicago

  • e[i]: the weight in the relationship between crime rate and percentage of nonewhite people in Chicago

  • percent_nonwhite[i]: the percentage of nonwhite people in zipcode i.

Specify the Model

This model resembles Model 3 but now our prediction of the crime rate accounts for both income levels and the percentage of nonwhites in the area.

Data

\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \alpha + u_{i} + \delta_{i} * time_{j} + \sigma_{i} * income_{i} + \rho_{i} * percentNonwhite_{i}\]

Priors

\[\alpha \sim N(0, 1000^2)\] \[\beta \sim N(0, 1000^2)\] \[u_{i} \sim N(0, prec.u)\] \[\delta_{i} \sim N(0, prec.delta)\] \[\sigma_{i} \sim N(0, prec.sigma)\] \[\rho_{i} \sim N(0, 0.01)\]

Hyperprior

\[prec.u \sim Gamma(0.5, 0.0005)\] \[prec.delta \sim Gamma(0.5, 0.0005)\] \[prec.sigma \sim Gamma(0.5, 0.0005)\]

#Specify the model
crime_model_4 <- "model{
    #Data
    for(i in 1:58) {
      for(j in 1:8){
        y[i,j] ~ dbin(p[i,j], n[i])
        logit(p[i,j]) = delta[i] * time[j] + sigma[i]*income[i] + rho[i]*percent_nonwhite[i]
      }
    }
    
    #Priors
    for(i in 1:58){
      delta[i] ~ dnorm(0, prec.delta)
      sigma[i] ~ dnorm(0, prec.sigma)
      rho[i] ~ dnorm(0, 0.01)
    }

    #Hyperpriors
    prec.delta ~ dgamma(0.5, 0.0005)
    prec.sigma ~ dgamma(0.5, 0.0005)
    
}"

#set up an algorithm to simulate the posterior by 
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_4 <- jags.model(textConnection(crime_model_4), data=list(y=Crime_count, n=Population.vector, time=time, income=Model_Data$per_capita_income, percent_nonwhite=100 - Model_Data$percent_white),
                    inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 464
##    Unobserved stochastic nodes: 176
##    Total graph size: 3144
## 
## Initializing model
#simulate a sample from the posterior 
#note that we specify both mu and tau variables
crime_sim_4 <- coda.samples(crime_jags_4, variable.names=c("p", "delta", "sigma", "rho"), n.iter=40000)

#store the samples in a data frame:    
crime_data_4 <- data.frame(crime_sim_4[[1]])

Running Mean Plots

running_mean_plot(crime_data_4$delta.4.)

running_mean_plot(crime_data_4$sigma.4.)

running_mean_plot(crime_data_4$rho.4.)

Crime Model 4 Visualizations

crime_viz_4_delta <- crime_data_4 %>%
  colMeans() %>%
  as.data.frame() %>%
  rownames_to_column(var="zipcode") %>%
  setNames(., c("zipcode", "delta")) %>%
  filter(grepl("delta", zipcode)) %>%
  mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
  left_join(zip_data, by="zipcode") %>%
  mutate(id=as.character(zip))
crime_model_4_delta <- left_join(zip_df, crime_viz_4_delta, by="id")
head(crime_model_4_delta)
##        long      lat order  hole piece    id   group zipcode       delta
## 1 -87.62271 41.88884     1 FALSE     1 60601 60601.1       1 -0.02008539
## 2 -87.62232 41.88879     2 FALSE     1 60601 60601.1       1 -0.02008539
## 3 -87.62120 41.88866     3 FALSE     1 60601 60601.1       1 -0.02008539
## 4 -87.62064 41.88860     4 FALSE     1 60601 60601.1       1 -0.02008539
## 5 -87.62058 41.88859     5 FALSE     1 60601 60601.1       1 -0.02008539
## 6 -87.62055 41.88859     6 FALSE     1 60601 60601.1       1 -0.02008539
##     zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
crime_viz_4_sigma <- crime_data_4 %>%
  colMeans() %>%
  as.data.frame() %>%
  rownames_to_column(var="zipcode") %>%
  setNames(., c("zipcode", "sigma")) %>%
  filter(grepl("sigma", zipcode)) %>%
  mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
  left_join(zip_data, by="zipcode") %>%
  mutate(id=as.character(zip))
crime_model_4_sigma <- left_join(zip_df, crime_viz_4_sigma, by="id")
head(crime_model_4_sigma)
##        long      lat order  hole piece    id   group zipcode        sigma
## 1 -87.62271 41.88884     1 FALSE     1 60601 60601.1       1 0.0001368507
## 2 -87.62232 41.88879     2 FALSE     1 60601 60601.1       1 0.0001368507
## 3 -87.62120 41.88866     3 FALSE     1 60601 60601.1       1 0.0001368507
## 4 -87.62064 41.88860     4 FALSE     1 60601 60601.1       1 0.0001368507
## 5 -87.62058 41.88859     5 FALSE     1 60601 60601.1       1 0.0001368507
## 6 -87.62055 41.88859     6 FALSE     1 60601 60601.1       1 0.0001368507
##     zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
crime_viz_4_rho <- crime_data_4 %>%
  colMeans() %>%
  as.data.frame() %>%
  rownames_to_column(var="zipcode") %>%
  setNames(., c("zipcode", "rho")) %>%
  filter(grepl("rho", zipcode)) %>%
  mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
  left_join(zip_data, by="zipcode") %>%
  mutate(id=as.character(zip))
crime_model_4_rho <- left_join(zip_df, crime_viz_4_rho, by="id")
head(crime_model_4_rho)
##        long      lat order  hole piece    id   group zipcode       rho
## 1 -87.62271 41.88884     1 FALSE     1 60601 60601.1       1 0.5684385
## 2 -87.62232 41.88879     2 FALSE     1 60601 60601.1       1 0.5684385
## 3 -87.62120 41.88866     3 FALSE     1 60601 60601.1       1 0.5684385
## 4 -87.62064 41.88860     4 FALSE     1 60601 60601.1       1 0.5684385
## 5 -87.62058 41.88859     5 FALSE     1 60601 60601.1       1 0.5684385
## 6 -87.62055 41.88859     6 FALSE     1 60601 60601.1       1 0.5684385
##     zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601

Map of Crime Differential Growth by Zipcode (Using Model 4)

model_4_map <- ggplot() 
model_4_map <- model_4_map +  geom_polygon(data = crime_model_4_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2) 
model_4_map <- model_4_map + coord_map() 
model_4_map <- model_4_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_4_map <- model_4_map + cleanup 
print(model_4_map)

Map of Income Differential Growth by Zipcode (Using Model 4)

model_4_map_sigma <- ggplot() 
model_4_map_sigma <- model_4_map_sigma +  geom_polygon(data = crime_model_4_sigma, aes(x=long, y=lat, group=group, fill=sigma), color = "black", size=0.2) 
model_4_map_sigma <- model_4_map_sigma + coord_map() 
model_4_map_sigma <- model_4_map_sigma + labs(x="", y="", title=paste("Income Differential Growth by Zipcode"), fill="total")
model_4_map_sigma <- model_4_map_sigma + cleanup
print(model_4_map_sigma)

Map of Percentage Nonwhite Differential Growth by Zipcode (Using Model 4)

model_4_map_rho <- ggplot() 
model_4_map_rho <- model_4_map_rho +  geom_polygon(data = crime_model_4_rho, aes(x=long, y=lat, group=group, fill=rho), color = "black", size=0.2) 
model_4_map_rho <- model_4_map_rho + coord_map() 
model_4_map_rho <- model_4_map_rho + labs(x="", y="", title=paste("Percentage Nonwhite Differential Growth by Zipcode"), fill="total")
model_4_map_rho <- model_4_map_rho + cleanup
print(model_4_map_rho)


Police Traffic Stops model

Cleaning the data

Data frame for the total population by race and zipcode

#Number of zipcodes: 62; Number of races: 7
Police <- read.csv("Police_stops_2016_2017_sample.csv")

Police.r <- Police %>%
  filter(!(RACE_CODE_CD=="I" | RACE_CODE_CD=="P" | RACE_CODE_CD=="WHI") ) %>%
  mutate(RACE_CODE_CD = as.character(RACE_CODE_CD)) %>%
  mutate(RACE_CODE_CD=replace(RACE_CODE_CD, RACE_CODE_CD=="WWH", "HISPANIC"))

vars = c("percent_white","percent_black","percent_asian","percent_hispanic")
Population.race <- Model_Data %>%
  mutate_at(vars, function(x){trunc(Model_Data$total_population*x/100)}) %>%
  select(zip, percent_asian,percent_black,percent_hispanic, percent_white) %>%
  unique() %>%
  remove_rownames %>%
  column_to_rownames(var='zip')
head(Population.race)
##       percent_asian percent_black percent_hispanic percent_white
## 60601          1751          1070              583          6131
## 60602           321            29              175           907
## 60603           264            35               44           538
## 60604            66             0               12           332
## 60605          4245          4494             1248         14483
## 60606           362            55              251          2091

Data frame for the number of arrests by zipcode and race

valid_zips <- unique(rownames(Population.race))
Arrests <- Police.r %>% 
  filter(ZIP_CD %in% valid_zips) %>%
  group_by(ZIP_CD, RACE_CODE_CD) %>%
  dplyr::summarise(total.arrest = as.integer(sum(ENFORCEMENT_ACTION_TAKEN_I))) %>%
  spread(key=ZIP_CD, value=total.arrest) %>%
  as.data.frame() %>%
  remove_rownames %>%
  column_to_rownames(var='RACE_CODE_CD')
Arrests[is.na(Arrests)] <- as.integer(0)
head(Arrests)
##          60601 60602 60603 60604 60605 60606 60607 60608 60609 60610 60611
## API          0     0     0     0     0     0     0     0     1     1     0
## BLK          7     2     2     3     5     2    13    13    87    41     3
## HISPANIC     3     0     1     0     1     0     1    39    43     1     0
## WHT          2     0     0     0     4     0     5     7     8     3     0
##          60612 60613 60614 60615 60616 60617 60618 60619 60620 60621 60622
## API          0     1     0     0     3     0     0     0     1     0     0
## BLK         76     9     7    39    16   100     3    88   103    97     7
## HISPANIC     6     6     3     0     1    19    28     0     2     0    10
## WHT          9     7     4     0     3     3    11     0     3     0     6
##          60623 60624 60625 60626 60628 60629 60630 60631 60632 60633 60634
## API          0     2     1     0     4     0     0     0     0     0     0
## BLK        147   352     6    31   121    40     1     1     5     0     3
## HISPANIC    27     9    21     7     0    27     3     1    40     1     8
## WHT          7    12     4     2     1     4     7     3     5     0     4
##          60636 60637 60638 60639 60640 60641 60642 60643 60644 60645 60646
## API          1     0     0     0     1     0     0     0     0     0     0
## BLK        128    92     1     9    29     4     1    25   211     4     0
## HISPANIC     2     1     1    25     7    13     3     0     3     7     0
## WHT          2     0     4     4     7     8     1     2     9     3     1
##          60647 60649 60651 60652 60653 60654 60655 60656 60657 60659 60660
## API          0     0     0     0     0     0     0     0     0     3     2
## BLK          4    63    96     7    42     6     0     0    10     1    11
## HISPANIC    23     1    29     5     0     0     1     1     4     4     1
## WHT          3     1     6     2     0     1     2     0     2     6     2
##          60661 60707 60827
## API          0     0     0
## BLK          1     2     7
## HISPANIC     0     3     0
## WHT          1     2     0

Data frame for the number of police stops by zipcode and race

Stops <- Police.r %>%
  filter(ZIP_CD %in% valid_zips) %>%
  group_by(ZIP_CD, RACE_CODE_CD) %>%
  dplyr::summarise(total.stop = n()) %>%
  spread(key=ZIP_CD, value=total.stop) %>%
  as.data.frame() %>%
  remove_rownames %>%
  column_to_rownames(var='RACE_CODE_CD')
Stops[is.na(Stops)] <- 0
head(Stops)
##          60601 60602 60603 60604 60605 60606 60607 60608 60609 60610 60611
## API          0     0     0     0     0     0     0     0     2     2     0
## BLK         15     4     6    10    16     6    36    50   291   121     9
## HISPANIC     4     0     1     1     4     0     2   192   187     7     3
## WHT          4     0     2     1     5     0    11    27    36    11     4
##          60612 60613 60614 60615 60616 60617 60618 60619 60620 60621 60622
## API          1     1     1     0     4     1     0     1     2     0     1
## BLK        245    37    24   154    63   351    15   296   339   339    31
## HISPANIC    23    21    10     2     8    80    60     2     5     2    37
## WHT         24    24    12     4     8    13    22     1     5     2    21
##          60623 60624 60625 60626 60628 60629 60630 60631 60632 60633 60634
## API          1     4     6     2     4     3     8     0     1     0     0
## BLK        426   899    16   110   414   184    17     3    19     1     9
## HISPANIC   134    43    50    44     8   130    29     2   165     2    32
## WHT         21    71    32    27     9    29    64     9    28     0    23
##          60636 60637 60638 60639 60640 60641 60642 60643 60644 60645 60646
## API          3     2     2     0     5     0     0     0     1     3     1
## BLK        389   299     9    37   111    13    13    86   626    38     1
## HISPANIC    12     1    12   107    23    51    10     1    23    17     3
## WHT          6     5    18    17    25    25     5     9    25    13     4
##          60647 60649 60651 60652 60653 60654 60655 60656 60657 60659 60660
## API          1     4     0     0     1     0     0     1     0     8     5
## BLK         19   237   320    32   203    18     2     0    46    24    32
## HISPANIC    87     3    95    11     2     0     4     1    16    31    14
## WHT          8     5    41    10     4     7     7     5    19    23     7
##          60661 60707 60827
## API          1     0     0
## BLK          3     3    31
## HISPANIC     0     7     1
## WHT          4     2     0

Preliminary visualizations

Stop Rate by Race

total_pop <- Model_Data %>% filter(Year==2010) %>% select(total_population)
population_by_ethnic <- Model_Data %>%
  filter(Year==2010) %>%
  mutate_at(c("percent_asian", "percent_black", "percent_hispanic", "percent_white"), function(x){x* total_pop$total_population/100}) %>%
  mutate(zip = as.numeric(zip)) %>%
  select(percent_asian, percent_black, percent_hispanic, percent_white, zip)
population_by_ethnic
##    percent_asian percent_black percent_hispanic percent_white   zip
## 1        1751.76       1070.52           583.92       6131.16 60601
## 2         321.86         29.26           175.56        907.06 60602
## 3         264.60         35.28            44.10        538.02 60603
## 4          66.40          0.00            12.45        332.00 60604
## 5        4245.24       4494.96          1248.60      14483.76 60605
## 6         362.57         55.78           251.01       2091.75 60606
## 7        4775.46       4021.44          2010.72      13321.02 60607
## 8        7879.40      12607.04         45700.52      11819.10 60608
## 9        2510.92      17576.44         30758.77      10671.41 60609
## 10       2711.80       6585.80          2324.40      26343.20 60610
## 11       4458.30       1188.88          1188.88      22291.50 60611
## 12       1442.08      21991.72          5047.28       7210.40 60612
## 13       2972.16       2972.16          4953.60      37647.36 60613
## 14       3418.15       2734.52          4785.41      56057.66 60614
## 15       2852.15      24447.00          2037.25      10186.25 60615
## 16      19301.49      11877.84          5444.01      11877.84 60616
## 17          0.00      45054.90         30870.95       5840.45 60617
## 18       3896.52       1948.26         45784.11      42861.72 60618
## 19          0.00      62974.34           649.22        649.22 60619
## 20          0.00      71055.41           732.53        732.53 60620
## 21          0.00      32263.68           336.08        336.08 60621
## 22       1628.16       4341.76         15196.16      32020.48 60622
## 23          0.00      28953.72         53649.54       1703.16 60623
## 24          0.00      37318.85           785.66        785.66 60624
## 25      11055.24       3158.64         30007.08      33165.72 60625
## 26       3574.27      12765.25         11744.03      21445.62 60626
## 27          0.00      67407.40          2151.30       1434.20 60628
## 28       1138.33      26181.59         76268.11      10244.97 60629
## 29       6367.46        578.86         15629.22      34152.74 60630
## 30        577.98        288.99          3467.88      24275.16 60631
## 31       4542.25       1816.90         75401.35       9084.50 60632
## 32        134.74       3368.50          5659.08       4311.68 60633
## 33       2908.72        727.18         25451.30      42903.62 60634
## 34          0.00      37522.92          1197.54        399.18 60636
## 35       1930.80      37167.90           965.40       7723.20 60637
## 36        560.33       1680.99         22973.53      30257.82 60638
## 37          0.00      13790.10         70789.18       6435.38 60639
## 38       8367.45      11585.70          8367.45      34113.45 60640
## 39       2829.20       1414.60         38194.20      26877.40 60641
## 40        915.60       2014.32          3845.52      10804.08 60642
## 41          0.00      37321.92          1555.08      11922.28 60643
## 42          0.00      44773.14          1428.93        952.62 60644
## 43       7036.95       7506.08         10320.86      20641.72 60645
## 44       2831.40        283.14          3397.68      21518.64 60646
## 45       2614.38       5228.76         45315.92      33115.48 60647
## 46          0.00      42717.36           908.88        908.88 60649
## 47          0.00      38524.50         20179.50       1834.50 60651
## 48          0.00      20408.16         14880.95       6802.72 60652
## 49        310.86      28288.26           621.72       1243.44 60653
## 50       1588.50        794.25          1429.65      11596.05 60654
## 51        281.73       1972.11          1690.38      24228.78 60655
## 52       2197.44        274.68          3296.16      21150.36 60656
## 53       4766.16       2042.64          4766.16      55151.28 60657
## 54      12690.24       1982.85          7931.40      15862.80 60659
## 55       5467.80       6309.00          7570.80      21450.60 60660
## 56       1635.92        446.16           446.16       4759.04 60661
## 57       1281.36       2989.84         12386.48      24772.96 60707
## 58          0.00      26263.24          1141.88        856.41 60827
print(colSums(population_by_ethnic))
##    percent_asian    percent_black percent_hispanic    percent_white 
##         152464.1         885195.5         785951.5         893236.1 
##              zip 
##        3516807.0
data_stops = data.frame(race=c("Asian Pacific Islander", "Black", "Hispanic", "White"), stop_rate=(rowSums(Stops)/colSums(population_by_ethnic %>% select(-zip))))
rownames(data_stops) <- c()
head(data_stops %>% arrange(desc(stop_rate)))
##                     race    stop_rate
## 1                  Black 0.0080750524
## 2               Hispanic 0.0023182092
## 3                  White 0.0009448790
## 4 Asian Pacific Islander 0.0005443903

Blacks and Hispanics are stopped at higher rate than Whites and Asians given their population.

Plot of Traffic Stops by Population

l <- population_by_ethnic %>%
  mutate(BLK=percent_black, HISPANIC=percent_hispanic, API=percent_asian, WHT=percent_white) %>%
  select(BLK, HISPANIC, API, WHT, zip) %>%
  gather(key=RACE_CODE_CD, value=population, - zip)

data <- Police.r %>%
  filter(ZIP_CD %in% valid_zips) %>%
  group_by(ZIP_CD, RACE_CODE_CD) %>%
  dplyr::summarise(total.stop = n()) %>%
  inner_join(l, by=c("ZIP_CD" = 'zip', 'RACE_CODE_CD')) %>%
  mutate(race=RACE_CODE_CD)

ggplot(data, aes(x=population, y=total.stop, colour=RACE_CODE_CD)) + geom_point() + geom_smooth() +
  labs(x="Population Count", y="Number of Traffic Stops", title="Traffic Stops By Population", subtitle="Visualization of number of traffic stops as the population of each ethnic group increases for each zipcode") + 
  guides(colour=guide_legend("Race")) +
  theme(plot.title = element_text(hjust=0.5)) 
## `geom_smooth()` using method = 'loess'

Plot of Traffic Stops by Percentage of Population

ggplot(data, aes(x=population, y=total.stop, colour=RACE_CODE_CD)) + geom_point() + geom_smooth() +
  labs(x="Percent of Population", y="Number of Traffic Stops", fill="Race", title="Traffic Stops By Percentage of Population") + 
  guides(colour=guide_legend("Race")) +
  theme(plot.title = element_text(hjust=0.5))
## `geom_smooth()` using method = 'loess'


Police Model 0

Variables Used

  • Data is only from January 2016 - February 2017 (13/12)

  • i: ethnic group

  • j: zipcode locations in Chicago

  • y[i, j]: observed number of arrests of ethnic group i in area j

  • n[i, j]: the number of traffic stops in area j for ethnic group i

  • \(\mu\): grand mean of arrests

  • \(\alpha\)[i]: the random effect for each ethnic group i

Data

\[y_{ij} \sim Pois(\frac{13}{12} n_{ij} \theta_{ij})\] \[log(\theta_{ij}) = \mu + \alpha_{i} \]

Priors

\[\mu \sim N(0, 0.000001)\] \[\alpha_{i} \sim N(0, 0.000001)\]

log(0.28603)
## [1] -1.251659
log(0.28603+0.05911) - log(0.28603)
## [1] 0.1878534
log(0.28603-0.05712) - log(0.28603)
## [1] -0.2227678
log(0.28603-0.12381) - log(0.28603)
## [1] -0.5671433
police_model_0 <- "model{
  #Data
  for (i in 1:4) {
    for(j in 1:58){
      y[i,j] ~ dpois((13/12) * n[i,j] * theta[i,j])
      log(theta[i,j]) <- mu + alpha[i]
    }
  }

  #Priors
  mu ~ dnorm(0, 1.0E-06)

  for (i in 1:4){
    alpha[i] ~ dnorm(0, 0.005)
  }
}"

#initial value set
inits = list(mu=-1.251659, alpha=c(0, 0.1878534, -0.2227678, -0.5671433))

police_jags_0 <- jags.model(textConnection(police_model_0), data=list(y=Arrests, n=Stops),
                    inits=inits)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 232
##    Unobserved stochastic nodes: 5
##    Total graph size: 597
## 
## Initializing model
#simulate a sample from the posterior 
#note that we specify both mu and tau variables
police_sim_0 <- coda.samples(police_jags_0, variable.names=c("mu", "alpha"), n.iter=40000)

#store the samples in a data frame:    
police_data_0 <- data.frame(police_sim_0[[1]])

Police Model 1

#estimate initial values
data <- Police.r %>% 
  filter(ZIP_CD %in% valid_zips) %>%
  group_by(ZIP_CD, RACE_CODE_CD) %>%
  dplyr::summarise(total.arrest = as.integer(sum(ENFORCEMENT_ACTION_TAKEN_I)), total.stop = n())
head(data)
## # A tibble: 6 x 4
## # Groups:   ZIP_CD [3]
##   ZIP_CD RACE_CODE_CD total.arrest total.stop
##    <int>        <chr>        <int>      <int>
## 1  60601          BLK            7         15
## 2  60601     HISPANIC            3          4
## 3  60601          WHT            2          4
## 4  60602          BLK            2          4
## 5  60603          BLK            2          6
## 6  60603     HISPANIC            1          1
freq_model <- lm(total.arrest ~ RACE_CODE_CD * total.stop, data = data)
summary(freq_model)
## 
## Call:
## lm(formula = total.arrest ~ RACE_CODE_CD * total.stop, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.097  -1.148  -0.147   1.273  46.686 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                     -0.08564    1.60786  -0.053  0.95758   
## RACE_CODE_CDBLK                 -4.88013    1.84214  -2.649  0.00875 **
## RACE_CODE_CDHISPANIC             0.54742    1.85234   0.296  0.76791   
## RACE_CODE_CDWHT                  1.12422    1.95807   0.574  0.56655   
## total.stop                       0.28603    0.49111   0.582  0.56098   
## RACE_CODE_CDBLK:total.stop       0.05911    0.49113   0.120  0.90433   
## RACE_CODE_CDHISPANIC:total.stop -0.05712    0.49135  -0.116  0.90757   
## RACE_CODE_CDWHT:total.stop      -0.12381    0.49393  -0.251  0.80235   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.55 on 189 degrees of freedom
## Multiple R-squared:  0.9789, Adjusted R-squared:  0.9781 
## F-statistic:  1250 on 7 and 189 DF,  p-value: < 2.2e-16

Variables Used

  • Data is only from January 2016 - February 2017 (13/12)

  • i: ethnic group

  • j: zipcode locations in Chicago

  • y[i, j]: observed number of arrests of ethnic group i in area j

  • n[i, j]: the number of traffic stops in area j for ethnic group i

  • \(\mu\): grand mean of arrests

  • \(\alpha\)[i]: the random effect for each ethnic group i

  • \(\beta\)[j]: the random effect for each zipcode location j

  • \(\epsilon\)[i, j]: the random error for each ethnic group i and location j

Data

\[y_{ij} \sim Pois(\frac{13}{12} n_{ij} \theta_{ij})\] \[log(\theta_{ij}) = \mu + \alpha_{i} + \beta_{j} + \epsilon_{ij}\]

Priors

\[\beta_{0} \sim N(0, 0.000001)\] \[\alpha_{i} \sim N(0, 0.000001)\] \[\beta_{j} \sim N(0, b^2)\] \[\epsilon_{ij} \sim N(0, e^2)\]

Hyperpriors

\[b \sim Gamma(0.5, 0.0005)\] \[e \sim Gamma(0.5, 0.0005)\]

police_model_1 <- "model{
  #Data
  for (i in 1:4) {
    for(j in 1:58){
      y[i,j] ~ dpois((13/12)*n[i,j]*theta[i,j])
      log(theta[i,j]) <- mu + alpha[i] + beta[j] + epsilon[i,j]
    }
  }

  #Priors

  mu ~ dnorm(0,1.0E-06)

  for (i in 1:4){
    alpha[i] ~ dnorm(0, 0.005)
  }

  for (j in 1:58){
    beta[j] ~ dnorm(0, 1/b^2)
  }

  for (i in 1:4){
    for (j in 1:58){
      epsilon[i,j] ~ dnorm(0, 1/e^2)
    }
  }

  b ~ dgamma(0.5, 0.0005)
  e ~ dgamma(0.5, 0.0005)
  
}"

#initial value set
inits = list(mu=-1.251659, alpha=c(0, 0.1878534, -0.2227678, -0.5671433), tau_beta=170, tau_epsilon=170, beta=rep(0.03, 58), epsilon=matrix(rep(0.01, 58*4), ncol=58))

police_jags_1 <- jags.model(textConnection(police_model_1), data=list(y=Arrests, n=Stops),
                    inits=inits)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 232
##    Unobserved stochastic nodes: 297
##    Total graph size: 1471
## 
## Initializing model
#simulate a sample from the posterior 
#note that we specify both mu and tau variables
police_sim_1 <- coda.samples(police_jags_1, variable.names=c("mu", "alpha", "beta", "sd_beta", "sd_epsilon", "epsilon"), n.iter=40000)

#store the samples in a data frame:    
police_data_1 <- data.frame(police_sim_1[[1]])

Running mean plot of beta.1., mu, and alpha.2.

running_mean_plot(police_data_1$beta.1.)

running_mean_plot(police_data_1$mu)

running_mean_plot(police_data_1$alpha.2.)

police_viz <- police_data_1 %>%
  colMeans() %>%
  as.data.frame() %>%
  rownames_to_column(var="zipcode") %>%
  setNames(., c("zipcode", "beta")) %>%
  filter(grepl("beta", zipcode)) %>%
  mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
  mutate(mu = mean(police_data_1$mu)) %>%
  mutate(rate = exp(mu + beta)) %>%
  left_join(zip_data, by="zipcode") %>%
  mutate(id=as.character(zip))
police_model_beta <- left_join(zip_df, police_viz, by="id")
head(police_model_beta)
##        long      lat order  hole piece    id   group zipcode       beta
## 1 -87.62271 41.88884     1 FALSE     1 60601 60601.1       1 0.05094085
## 2 -87.62232 41.88879     2 FALSE     1 60601 60601.1       1 0.05094085
## 3 -87.62120 41.88866     3 FALSE     1 60601 60601.1       1 0.05094085
## 4 -87.62064 41.88860     4 FALSE     1 60601 60601.1       1 0.05094085
## 5 -87.62058 41.88859     5 FALSE     1 60601 60601.1       1 0.05094085
## 6 -87.62055 41.88859     6 FALSE     1 60601 60601.1       1 0.05094085
##           mu      rate   zip
## 1 -0.2179287 0.8462099 60601
## 2 -0.2179287 0.8462099 60601
## 3 -0.2179287 0.8462099 60601
## 4 -0.2179287 0.8462099 60601
## 5 -0.2179287 0.8462099 60601
## 6 -0.2179287 0.8462099 60601

Map of rate at which traffic stops result in arrests by zipcode (Using Model 3)

model_1_map_beta <- ggplot() 
model_1_map_beta <- model_1_map_beta +  geom_polygon(data = police_model_beta, aes(x=long, y=lat, group=group, fill=rate), color = "black", size=0.2) 
model_1_map_beta <- model_1_map_beta + coord_map() 
model_1_map_beta <- model_1_map_beta + labs(x="", y="", title=paste("Rate at Which Traffic Stops Result in Arrests by Zipcode"), fill="total")
model_1_map_beta <- model_1_map_beta + cleanup
print(model_1_map_beta)

Rate of traffic stops by ethnic group

police_viz1 <- police_data_1 %>%
  colMeans() %>%
  as.data.frame() %>%
  rownames_to_column(var="zipcode") %>%
  setNames(., c("zipcode", "alpha")) %>%
  filter(grepl("alpha", zipcode)) %>%
  mutate(mu = mean(police_data_1$mu)) %>%
  mutate(rate = exp(mu + alpha))
head(police_viz1)
##    zipcode     alpha         mu      rate
## 1 alpha.1. -1.255402 -0.2179287 0.2291610
## 2 alpha.2. -1.097255 -0.2179287 0.2684249
## 3 alpha.3. -1.279364 -0.2179287 0.2237350
## 4 alpha.4. -1.344437 -0.2179287 0.2096395